% =========================================================================
% lmj_parperturb
% 1. Load in
% a) one  parameter mode and perturb one coefficient at a time 
%    according to partvec. True value will be added automatically 
% b) a table of modes 
% Also load the dataset 
% =========================================================================
clear all;close all;clc;cucd=pwd; 

% =========================================================================
%% user defined inputs 
% 1. What to run? 
% ==1  perturb one element at a time
% ==2  load a matrix of coefficients 
switch compare_type
    case 1
        parf_path = 'O:\PROJ_LIB\LJM\exog\udnosepmebnd\lindtr\1trend z tight u addL2Y fmincon'; 
        parf_name = 'top mode';
        partab_name = 'Sheet1'; 
        par_pert = 'alpha'; 
        pertvec = [0.15 0.25]; 
        savefolder = ['par_pert\', par_pert] ; 
    case 2
        parf_path = 'O:\PROJ_LIB\LJM\exog\udnosepmebnd\lindtr\1trend z tight u addL2Y fmincon'; 
        parf_name = 'tight';
        partab_name = 'modes'; 
        savefolder = ['compare sepnobs', strdate]; 
end
% 2. Data settings 
dataf_path = 'O:\PROJ_LIB\LJM\exog\data\April 2010';
dataf_name = 'data_detrendednosep'; 



funcmod = @modudnosepmebndlindtr;

    


ana_in.state_sel ={'yobs', 'cobs', 'iobs','uobs', 'fobs', 'lam','vobs','h','ls','L'}; % state selected
ana_in.flag_adds = [   1,      1,      1,   0,     0,       0,    0,     0    0   0]; % transform 2 levels
ana_in.specrep   ={'yobs lev','cobs lev','iobs lev','uobs','fobs','lam','vobs','h','ls','L'}; 
ana_in.lowper    =6; 
ana_in.highper   =32; 
ana_in.npoints   =2000; 

ssreport = {'u', 'f', 'L2Y100','L2Yout','L2Ylab','lshare', 'V', 'C2Y', 'I2Y', 'p_omegalow', 'p_omegahigh'}; 


%% Additional inputs
% =========================================================================
% 1. OPTIONS for  FZERO or FSOLVE solution of the model 
% TolFun will be taken as the tolerance criteria 
% =========================================================================
solveopt=optimset('MaxFunEvals',20000,'Display','off',...
    'MaxIter',5000,'TolFun',1e-8);

% =========================================================================
% 2. ADDSOL 
% Structure with additional inputs for the solution of the model 
% Only one required is .JFRVEC which is the initial guess grid for a solution 
% =========================================================================
addsol.jfrvec=[(0.00001:0.01:0.4)';(0.405:0.005:0.9)';(0.901:0.01:0.999)']; 

% =========================================================================
% 3. Data Setting
% =========================================================================
dataopt.start=1967.25;
dataopt.starest=1967.25;
dataopt.finish=2009.75; 
dataopt.finishest=2009.75; 
dataopt.demean   =0; 
dataopt.demean_type=[]; 
onscreen = 0; 

flag_plotst = 1;  
flag_plotss = 1; 

%% Load in Parameter
cd(parf_path); 
[parmode, parnames] = xlsread(parf_name, partab_name); 
cd(cucd); 

numpar_loaded = size(parmode, 2);
if compare_type == 1 && numpar_loaded ~= 1; 
    error('only one parameter vector allowed under PERTURB mode');
end 

if compare_type == 2 && numpar_loaded == 1; 
    error('at leaset two parameter vector required under COMPARE mode'); 
end 

printcell([parnames num2cprec(parmode)]); 
quer('c'); 

% check par yields good results
for jj = 1:numpar_loaded
    [GG, RR, CONS, eu, SDX, ZZ, initss, ssvec, flag, ssnames, stanames ,shonames]...
        =feval(funcmod,parmode(:, jj),solveopt,addsol);
    if ~isequal(eu, [1; 1]);
        error('loaded parameter does not solve the model ')
    end
    
    
    znames=cell(size(ZZ,1),1);
    for ii=1:length(znames)
        temp      =find(ZZ(ii,:)==1 ) ;
        if length(temp) > 1
            error('Need 1 entry only for each row of ZZ')
        end
        znames(ii)=stanames( temp );
    end
end


%% Load in Data
cd(dataf_path); 
[Y,Ynames,sample,trainvec]=load_dset(dataf_name,dataf_path,dataopt,dataopt.demean_type,onscreen);
cd(cucd); 

%% solve for the model
%  Sorted in ascending order 

switch compare_type
    case 1
        parpos_pert = cellposition(par_pert, parnames);
        parpert_true = parmode(parpos_pert);
        pos_dupvar = find(abs(parpert_true - pertvec) < 1e-5);
        if ~isempty(pos_dupvar)
            pertvec(pos_dupvar) = nan;
            pertvec = pertvec(~isnan(pertvec));
        end
        pertvec=[parpert_true, pertvec];
        numpert = length(pertvec);
        pertmat = repmat(parmode, 1, numpert);
        pertmat(parpos_pert,:) = pertvec;
        % Create top cell
        modenames = cell(1, numpert);
        for ii = 1:numpert;
            if abs( pertvec(ii) - parpert_true ) < eps
                modenames{ii} = [par_pert, ' = ', num2str(parpert_true), ' (est)'];
            else
                modenames{ii} = [par_pert, ' = ', num2str(pertvec(ii))];
            end
        end
    case 2
        pertmat = parmode;
        numpert = numpar_loaded;
        modenames = cell(1, numpert);
        for ii = 1:numpert;
            modenames{ii} = ['mode ', num2str(ii)];
        end
end


%% 
ana_in.flag_addz = zeros(1,size(Y, 2));
ana_in.flag_block = 1; 
ana_in.stanames = stanames;
ana_in.shonames = shonames;
% ana_in.savepath = parf_path;
ana_in.dates = sample; 
ana_in.parnames = parnames; 
ana_in.vecsim = [100, size(Y,1)]; 
ana_in.nrep = 2000; 
ana_in.nper  =4;
ana_in.flag_acov = 0; 
ana_in.flag_unit =0;
ana_in.nirf      =16;
ana_in.nvdec     =40;
ana_in.flag_acdecomp=0;
ana_in.flag_medbonly = 1;
ana_in.savefolder = savefolder; 
ana_in.modelnames = modenames; 

ana_in.flag_plotst = flag_plotst; 
if compare_type == 2; 
    flag_plotss = 0; 
end 
ana_in.flag_plotss = flag_plotss; 
ana_in.ssreport = ssreport; 
ana_in.ssnames = ssnames; 

if compare_type == 1
    ana_in.pospert = parpos_pert; 
end

save_path = cr_dir(parf_path,savefolder); 
ana_in.savepath = save_path; 
switch compare_type
    case 1
        [out,outcell,hand_vec]=dsge_spectrum(GG,RR,ZZ,SDX,ana_in.specrep,stanames,shonames,znames,...
            ana_in.lowper,ana_in.highper,ana_in.npoints,save_path);
    
        [vdcom,vcvth,sdevf,ACzero,ACone,vcshocks,auone,vdech,vdech_st,vdcom_st]=vardecom_mode(GG,RR,SDX,ZZ,50);
        
        cd(save_path)
        xlswrite('variance decomposition tab',crtcell( vdcom, znames , shonames ,'Stationary Variance' ),'stationary' );
        cd(cucd);
    case 2
        save_path = cr_dir(parf_path,savefolder);
        outpath   = save_path; 
        
        ssvec_mat = zeros(length(ssnames), numpar_loaded); 
        
        for ii = 1:numpar_loaded
            [GG, RR, CONS, eu, SDX, ZZ, initss, ssvec, flag, ssnames, stanames ,shonames]...
                            =feval(funcmod,parmode(:, ii),solveopt,addsol);
            ssvec_mat(:, ii) = ssvec(:); 
            figname = ['mode', num2str(ii), '_specdecom'];             
            [out,spec_outcell,hand_vec]=dsge_spectrum(GG,RR,ZZ,SDX,ana_in.specrep,stanames,shonames,znames,...
            ana_in.lowper,ana_in.highper,ana_in.npoints,save_path, figname);
            spec_outcell = merge_cells({['Mode ', num2str(ii)]}, spec_outcell); 
        
            [vdcom,vcvth,sdevf,ACzero,ACone,vcshocks,auone,vdech,vdech_st,vdcom_st]=vardecom_mode(GG,RR,SDX,ZZ,50);
            var_outcell = crtcell( vdcom, znames , shonames ,['Mode ', num2str(ii),': Stationary Variance']); 
        
            if ii == 1;
                spec_outcell_full = spec_outcell; 
                var_outcell_full = var_outcell; 
            else
                spec_outcell_full = merge_cells( spec_outcell_full, spec_outcell, 2);
                var_outcell_full  = merge_cells( var_outcell_full,   var_outcell, 2); 
            end 
        end 
        
%         load workspace root spec subf loadpath outpath
        temp      = extrsubf(parf_path, cucd ); 
        tempos    = findstr(temp,'\'); 
        
        root=temp(tempos(1):tempos(2)-1); 
        spec=temp(tempos(2):tempos(3)-1); 
        subf=temp(tempos(3):end        ); 
        clear temp tempos; 
        
        outfolder = extrsubf(outpath, parf_path); 
        topcell = {'root', root; 
                   'spec', spec; 
                   'subf', subf; 
                   'outfolder', outfolder
                   'data', dataf_name};
               
        spec_outcell_full = merge_cells(topcell, spec_outcell_full, 2); 
        var_outcell_full = merge_cells(topcell, var_outcell_full, 2); 
        
        par_cell = [parnames, num2cprec(parmode, 10)]; 
        ss_cell = [ssnames, num2cprec(ssvec_mat)]; 
        mode_cell = merge_cells(par_cell, ss_cell, 2); 
        mode_cell = merge_cells(topcell, mode_cell, 2); 
        cd(save_path)
        xlswrite('compare',spec_outcell_full,'spectrum');
        xlswrite('compare', var_outcell_full,'stationary');
        xlswrite('compare', mode_cell,'modes');
        cd(cucd)
end 
close all
%ana_out = lmj_modeanalysis(funcmod, pertmat, Y, ana_in, solveopt, addsol); 


%% mode_analysis COMPUTE
ana_out = lmj_modeanalysis_compute(funcmod, pertmat, Y, ana_in, solveopt, addsol); 

% 1) plot volatilities
in_stdhist.modenames = ana_out.modelnames; 
in_stdhist.varnames  = znames; 
in_stdhist.savepath  = ana_in.savepath;  
in_stdhist.savename  = 'Z STD'; 
plot_stdhist(ana_out.stdsimz_mat, Y, in_stdhist)

% 2) plot irfs

% (2.1) plot all observables
in_irf.modenames = ana_out.modelnames; 
in_irf.varnames = znames; 
in_irf.shonames = shonames; 
in_irf.savepath = save_path; 
in_irf.savename = 'Z IRF'; 
in_irf.addlev = ana_in.flag_addz;
in_irf.flag_irfbysho = 1; 
plotmultiple_irfs(ana_out.irfz_mat, in_irf); 


% (2.2) plotting selected states. 
in_irf.modenames = ana_out.modelnames; 
in_irf.varnames =  ana_in.state_sel; 
in_irf.shonames = shonames; 
in_irf.savepath = save_path; 
in_irf.savename = 'Sel States IRF'; 
in_irf.addlev = ana_in.flag_adds; 
in_irf.flag_irfbysho = 1; 
plotmultiple_irfs(ana_out.irfs_mat, in_irf); 


% 3) plot Variance Decompositions.
    % 3.1) of the observables
in_vdec.modenames = ana_out.modelnames;
in_vdec.varnames =  znames;
in_vdec.shonames = shonames;
in_vdec.savepath = save_path;
in_vdec.savename = 'Z VDEC';
plot_vdecs(ana_out.vdechz_mat, in_vdec);
    
    
    % 3.2) of the variables
in_vdec.modenames = ana_out.modelnames;
in_vdec.varnames =  ana_in.state_sel;
in_vdec.shonames = shonames;
in_vdec.savepath = save_path;
in_vdec.savename = 'Sel State VDEC';
plot_vdecs(ana_out.vdechs_mat, in_vdec);


% 4) AUTOCORRELATION
    % 4.1) of the obervables; 
in_ac.modenames = ['data';  ana_out.modelnames(:)]; 
in_ac.varnames = znames;
in_ac.savepath = save_path; 
in_ac.savename = 'Obs AC Median'; 
    
ac_data=accov(Y,(0:ana_in.nper));
clear med_acsimz_mat;
for ii = 1:ana_out.num_valid + 1;
    if ii == 1;
        med_acsimz_mat(:,:,:,ii)=ac_data;
    else
        med_acsimz_mat(:,:,:,ii)=median(ana_out.acsimz_mat(:, :, :, :, ii-1),4);
    end
end
plotmultiple_acorr(med_acsimz_mat, in_ac);



% 4.2) of the selected states
in_ac.modenames = ana_out.modelnames; 
in_ac.varnames  = ana_in.state_sel; 
in_ac.savepath = save_path; 
in_ac.savename = 'Sel States AC Median'; 

clear med_acsims_mat
for ii = 1:ana_out.num_valid;
    med_acsims_mat(:,:,:,ii)=median(ana_out.acsims_mat(:, :, :, :, ii),4);
end
plotmultiple_acorr(med_acsims_mat, in_ac);


% 4.3) full cross-correlogram vs. data

in_acfull.modenames =  [ana_out.modelnames(:); 'data']; 
in_acfull.varnames  = znames; 
in_acfull.savepath  = save_path; 
in_acfull.savename = 'Obs AC full'; 
plotmultiple_acall(ana_out.acsimz_mat, ac_data, in_acfull); 


% 5, selected steady state values (if perturbing one parameter)
if compare_type == 1
    in_ssplot.modenames = ana_out.modelnames;
    in_ssplot.savepath  = save_path;
    in_ssplot.savename  = 'Sel Steady States';
    in_ssplot.ssreport  = ana_in.ssreport;
    in_ssplot.ssnames   = ssnames;
    in_ssplot.parnames  = parnames;
    in_ssplot.pertvar = par_pert;
    
    plotperturb_lmj(pertmat, funcmod, in_ssplot, solveopt, addsol);
end
